This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
#library(spdep)
#library(maptools)
#library(rgeos)
#library(GISTools)
library(sf)
## Warning: package 'sf' was built under R version 4.2.2
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
#library(sp)
#library(tmap)
#library(spatialreg)
#library(rgdal)
library(tidyverse)
## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.7 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(biscale)
## Warning: package 'biscale' was built under R version 4.2.2
#library(cowplot)
library(mapview)
## Warning: package 'mapview' was built under R version 4.2.2
setwd("C://Users//Kaifs//OneDrive//Documents//DHS individual//shpfile")
BJ <- st_read(dsn="BJGE71FL/BJGE71FL.shp") %>%
dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `BJGE71FL' from data source
## `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\BJGE71FL\BJGE71FL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 555 features and 20 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 5.684342e-14 ymin: 5.684342e-14 xmax: 3.611158 ymax: 12.34528
## Geodetic CRS: WGS 84
BU <- st_read(dsn="BUGE71FL/BUGE71FL.shp") %>%
dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `BUGE71FL' from data source
## `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\BUGE71FL\BUGE71FL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 554 features and 20 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 0 ymin: -4.432215 xmax: 30.784 ymax: 0
## Geodetic CRS: WGS 84
CM <- st_read(dsn="CMGE71FL/CMGE71FL.shp") %>%
dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `CMGE71FL' from data source
## `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\CMGE71FL\CMGE71FL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 430 features and 20 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 9.200917 ymin: 2.028929 xmax: 16.06679 ymax: 12.76956
## Geodetic CRS: WGS 84
GM <- st_read(dsn="GMGE81FL/GMGE81FL.shp") %>%
dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `GMGE81FL' from data source
## `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\GMGE81FL\GMGE81FL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 280 features and 20 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -16.80157 ymin: 13.14906 xmax: -13.84872 ymax: 13.80121
## Geodetic CRS: WGS 84
GN <- st_read(dsn="GNGE71FL/GNGE71FL.shp") %>%
dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `GNGE71FL' from data source
## `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\GNGE71FL\GNGE71FL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 401 features and 20 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -14.79391 ymin: 7.408167 xmax: -7.856225 ymax: 12.51171
## Geodetic CRS: WGS 84
NG <- st_read(dsn="NGGE7BFL/NGGE7BFL.shp") %>%
dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `NGGE7BFL' from data source
## `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\NGGE7BFL\NGGE7BFL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1389 features and 20 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 0 ymin: 0 xmax: 13.66089 ymax: 13.76644
## Geodetic CRS: WGS 84
LB <- st_read(dsn="LBGE7AFL/LBGE7AFL.shp") %>%
dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `LBGE7AFL' from data source
## `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\LBGE7AFL\LBGE7AFL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 325 features and 20 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -11.41053 ymin: 0 xmax: 0 ymax: 8.419803
## Geodetic CRS: WGS 84
ML <- st_read(dsn="MLGE7AFL/MLGE7AFL.shp") %>%
dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `MLGE7AFL' from data source
## `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\MLGE7AFL\MLGE7AFL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 345 features and 20 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -11.92693 ymin: 5.684342e-14 xmax: 1.42239 ymax: 18.45378
## Geodetic CRS: WGS 84
MR <- st_read(dsn="MRGE71FL/MRGE71FL.shp") %>%
dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `MRGE71FL' from data source
## `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\MRGE71FL\MRGE71FL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1200 features and 20 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -17.06289 ymin: 0 xmax: 0 ymax: 25.42278
## Geodetic CRS: WGS 84
RW <- st_read(dsn="RWGE81FL/RWGE81FL.shp") %>%
dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `RWGE81FL' from data source
## `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\RWGE81FL\RWGE81FL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 500 features and 20 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 28.89335 ymin: -2.8127 xmax: 30.82515 ymax: -1.078855
## Geodetic CRS: WGS 84
SL <- st_read(dsn="SLGE7AFL/SLGE7AFL.shp") %>%
dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `SLGE7AFL' from data source
## `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\SLGE7AFL\SLGE7AFL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 576 features and 20 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -13.28709 ymin: 0 xmax: 0 ymax: 9.976958
## Geodetic CRS: WGS 84
ZM <- st_read(dsn="ZMGE71FL/ZMGE71FL.shp") %>%
dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `ZMGE71FL' from data source
## `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\ZMGE71FL\ZMGE71FL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 545 features and 20 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -0.0048541 ymin: -17.94572 xmax: 33.60199 ymax: 0.00149906
## Geodetic CRS: WGS 84
map <- rbind(BJ,BU,CM,GM,GN,NG,LB,ML,MR,RW,SL,ZM)
#label var cum_pre_hdevi "cumulative abnormally high precipitation (m)"
#label var cum_pre_hdevi_E "cumulative abnormally high precipitation: early childhood(m)"
#label var cum_pre_hdevi_L "cumulative abnormally high precipitation: mid/late childhood (m)"
#label var cum_pre_ldevi "cumulative abnormally low precipitation (m)"#
#label var cum_pre_ldevi_E "cumulative abnormally low precipitation: early childhood(m)"
#label var cum_pre_ldevi_L "cumulative abnormally low precipitation: mid/late (m)"
#//temperature
#label var cum_mtemp_hdevi "cumulative abnormally high temperature (m)"
#label var cum_mtemp_hdevi_E "cumulative abnormally high temperature: early childhood(m)"
#label var cum_mtemp_hdevi_L "cumulative abnormally high temperature: mid/late childhood(m)"
#label var cum_mtemp_ldevi "cumulative abnormally low precipitation (m)"
#label var cum_mtemp_ldevi_E "cumulative abnormally low precipitation: early childhood (m)"
#label var cum_mtemp_ldevi_L "cumulative abnormally low precipitation: mid/late childhood (m)"
env_africa<-read.csv("TenAfricacountries.csv")
env_africa$DHSCC <- env_africa$dhscc
env_africa$DHSCLUST <- env_africa$dhsclust
env_africa<-env_africa %>%
group_by(dhscc,dhsclust,dhsyear,DHSCLUST,DHSCC,all_population_count_2015,global_human_footprint,nightlights_composite) %>%
summarise(cum_pre_hdevi_m = mean(cum_pre_hdevi),
cum_pre_hdevi_E_m = mean(cum_pre_hdevi_E),
cum_pre_hdevi_L_m = mean(cum_pre_hdevi_L),
cum_pre_ldevi_m = mean(cum_pre_ldevi),
cum_pre_ldevi_E_m = mean(cum_pre_ldevi_E),
cum_pre_ldevi_L_m = mean(cum_pre_ldevi_L),
cum_mtemp_hdevi_m = mean(cum_mtemp_hdevi),
cum_mtemp_hdevi_E_m = mean(cum_mtemp_hdevi_E),
cum_mtemp_hdevi_L_m = mean(cum_mtemp_hdevi_L),
cum_mtemp_ldevi_m = mean(cum_mtemp_ldevi),
cum_mtemp_ldevi_E_m = mean(cum_mtemp_ldevi_E),
cum_mtemp_ldevi_L_m = mean(cum_mtemp_ldevi_L),
average_prim = mean(pri_com1)) %>%
left_join(map, by=c("DHSCLUST","DHSCC"))
## `summarise()` has grouped output by 'dhscc', 'dhsclust', 'dhsyear', 'DHSCLUST',
## 'DHSCC', 'all_population_count_2015', 'global_human_footprint'. You can
## override using the `.groups` argument.
# right now try with point in time variable
#BJ_map <- unique( env_BJ[ ,c("DHSCLUST","all_population_count_2015","aridity_2015","average_prim")]) %>%
# left_join(BJ, by="DHSCLUST")
#CM_map <- unique( env_CM[ ,c("DHSCLUST","all_population_count_2015","aridity_2015","average_prim")]) %>%
# left_join(CM, by="DHSCLUST")
#GM_map <- unique( env_GM[ ,c("DHSCLUST","all_population_count_2015","aridity_2015","average_prim")]) %>%
# left_join(GM, by="DHSCLUST")
#NG_map <- unique( env_NG[ ,c("DHSCLUST","all_population_count_2015","aridity_2015","average_prim")]) %>%
# left_join(NG, by="DHSCLUST")
#ML_map <- unique( env_ML[ ,c("DHSCLUST","all_population_count_2015","aridity_2015","average_prim")]) %>%
# left_join(ML, by="DHSCLUST")
custom_pal3 <- c(
"1-1" = "#d3d3d3", # low x, low y
"2-1" = "#ba8890",
"3-1" = "#9e3547", # high x, low y
"1-2" = "#8aa6c2",
"2-2" = "#7a6b84", # medium x, medium y
"3-2" = "#682a41",
"1-3" = "#4279b0", # low x, high y
"2-3" = "#3a4e78",
"3-3" = "#311e3b" # high x, high y
)
legend<-bi_legend(pal = custom_pal3,
dim = 3,
xlab = "Higher Climate Hazards ",
ylab = "Higher Dropout ",
size = 12)
bi_class_pre_h <- c("1 - 1","2 - 1","3 - 1","1 - 2","2 - 2","3 - 2","1 - 3","2 - 3","3 - 3")
bi_class_pre_l <- c("1 - 1","2 - 1","3 - 1","1 - 2","2 - 2","3 - 2","1 - 3","2 - 3","3 - 3")
bi_class_tem_h <- c("1 - 1","2 - 1","3 - 1","1 - 2","2 - 2","3 - 2","1 - 3","2 - 3","3 - 3")
bi_class_tem_l <- c("1 - 1","2 - 1","3 - 1","1 - 2","2 - 2","3 - 2","1 - 3","2 - 3","3 - 3")
bivariate_color_scale <- data.frame(bi_class_pre_h,bi_class_pre_l,bi_class_tem_h,bi_class_tem_l,custom_pal3)
map_africa_pre_h <- env_africa %>%
group_by(dhscc) %>%
mutate(pre_h_class = as.numeric(cut(cum_pre_hdevi_m, 3,na.rm=T)),
dropout_class = case_when(average_prim<0.45 ~ 3,
average_prim<0.9 & average_prim>0.45 ~ 2,
average_prim >=0.9 ~ 1),
bi_class_pre_h = case_when(pre_h_class == 1 & dropout_class == 1 ~ "1 - 1",
pre_h_class == 2 & dropout_class == 1 ~ "1 - 2",
pre_h_class == 3 & dropout_class == 1 ~ "1 - 3",
pre_h_class == 1 & dropout_class == 2 ~ "2 - 1",
pre_h_class == 2 & dropout_class == 2 ~ "2 - 2",
pre_h_class == 3 & dropout_class == 2 ~ "2 - 3",
pre_h_class == 1 & dropout_class == 3 ~ "3 - 1",
pre_h_class == 2 & dropout_class == 3 ~ "3 - 2",
pre_h_class == 3 & dropout_class == 3 ~ "3 - 3")) %>%
ungroup() %>%
left_join(bivariate_color_scale, by = "bi_class_pre_h") %>%
filter(all_population_count_2015>0)
map_africa_pre_l <- env_africa %>%
group_by(dhscc) %>%
mutate(pre_h_class = as.numeric(cut(cum_pre_hdevi_m, 3,na.rm=T)),
pre_l_class = as.numeric(cut(cum_pre_ldevi_m, 3,na.rm=T)),
tem_h_class = as.numeric(cut(cum_mtemp_hdevi_m, 3,na.rm=T)),
tem_l_class = as.numeric(cut(cum_mtemp_ldevi_m, 3,na.rm=T)),
dropout_class = case_when(average_prim<0.45 ~ 3,
average_prim<0.9 & average_prim>0.45 ~ 2,
average_prim >=0.9 ~ 1),
bi_class_pre_l = case_when(pre_l_class == 1 & dropout_class == 1 ~ "1 - 1",
pre_l_class == 2 & dropout_class == 1 ~ "1 - 2",
pre_l_class == 3 & dropout_class == 1 ~ "1 - 3",
pre_l_class == 1 & dropout_class == 2 ~ "2 - 1",
pre_l_class == 2 & dropout_class == 2 ~ "2 - 2",
pre_l_class == 3 & dropout_class == 2 ~ "2 - 3",
pre_l_class == 1 & dropout_class == 3 ~ "3 - 1",
pre_l_class == 2 & dropout_class == 3 ~ "3 - 2",
pre_l_class == 3 & dropout_class == 3 ~ "3 - 3")) %>%
ungroup() %>%
left_join(bivariate_color_scale, by = "bi_class_pre_l") %>%
filter(all_population_count_2015>0)
map_africa_tem_h <- env_africa %>%
group_by(dhscc) %>%
mutate(tem_h_class = as.numeric(cut(cum_mtemp_hdevi_m, 3,na.rm=T)),
dropout_class = case_when(average_prim<0.45 ~ 3,
average_prim<0.9 & average_prim>0.45 ~ 2,
average_prim >=0.9 ~ 1),
bi_class_tem_h = case_when(tem_h_class == 1 & dropout_class == 1 ~ "1 - 1",
tem_h_class == 2 & dropout_class == 1 ~ "1 - 2",
tem_h_class == 3 & dropout_class == 1 ~ "1 - 3",
tem_h_class == 1 & dropout_class == 2 ~ "2 - 1",
tem_h_class == 2 & dropout_class == 2 ~ "2 - 2",
tem_h_class == 3 & dropout_class == 2 ~ "2 - 3",
tem_h_class == 1 & dropout_class == 3 ~ "3 - 1",
tem_h_class == 2 & dropout_class == 3 ~ "3 - 2",
tem_h_class == 3 & dropout_class == 3 ~ "3 - 3")) %>%
ungroup() %>%
left_join(bivariate_color_scale, by = "bi_class_tem_h") %>%
filter(all_population_count_2015>0)
map_africa_tem_l <- env_africa %>%
group_by(dhscc) %>%
mutate(tem_l_class = as.numeric(cut(cum_mtemp_ldevi_m, 3,na.rm=T)),
dropout_class = case_when(average_prim<0.45 ~ 3,
average_prim<0.9 & average_prim>0.45 ~ 2,
average_prim >=0.9 ~ 1),
bi_class_tem_l = case_when(tem_l_class == 1 & dropout_class == 1 ~ "1 - 1",
tem_l_class == 2 & dropout_class == 1 ~ "1 - 2",
tem_l_class == 3 & dropout_class == 1 ~ "1 - 3",
tem_l_class == 1 & dropout_class == 2 ~ "2 - 1",
tem_l_class == 2 & dropout_class == 2 ~ "2 - 2",
tem_l_class == 3 & dropout_class == 2 ~ "2 - 3",
tem_l_class == 1 & dropout_class == 3 ~ "3 - 1",
tem_l_class == 2 & dropout_class == 3 ~ "3 - 2",
tem_l_class == 3 & dropout_class == 3 ~ "3 - 3")) %>%
ungroup() %>%
left_join(bivariate_color_scale, by = "bi_class_tem_l") %>%
filter(all_population_count_2015>0)
mapview(na.omit(map_africa_pre_h),
zcol = "bi_class_pre_h",
xcol = "LONGNUM",
ycol = "LATNUM",
legend = F,
crs = 4269,
grid = FALSE,
#cex = 3,
cex="all_population_count_2015",
#alpha.regions = 0.2,
alpha = 0.5,
col.regions=list("#d3d3d3","#8aa6c2","#4279b0", "#ba8890","#7a6b84","#3a4e78","#9e3547","#682a41","#311e3b"))